Liquid-gas coexistence and critical point shifts in size-disperse fluids 
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Specialized Monte Carlo simulations and the moment free energy (MFE) method are employed 
to study liquid-gas phase equilibria in size-disperse fluids. The investigation is made subject to 
the constraint of fixed polydispersity, i.e. the form of the 'parent' density distribution fP{o) of 
the particle diameters cr, is prescribed. This is the experimentally realistic scenario for e.g. colloidal 
dispersions. The simulations are used to obtain the cloud and shadow curve properties of a Lennard- 
Jones fluid having diameters distributed according to a Schulz form with a large [5 ~ 40%) degree 
of polydispersity. Good qualitative accord is found with the results from a MFE method study 
of a corresponding van der Waals model that incorporates size-dispersity both in the hard core 
reference and the attractive parts of the free energy. The results show that polydispersity engenders 
considerable broadening of the coexistence region between the cloud curves. The principal effect 
of fractionation in this region is a common overall scaling of the particle sizes and typical inter- 
particle distances, and we discuss why this effect is rather specific to systems with Schulz diameter 
distributions. Next, by studying a family of such systems with distributions of various widths, we 
estimate the dependence of the critical point parameters on 5. In contrast to a previous theoretical 
prediction, size-dispersity is found to raise the critical temperature above its monodisperse value. 
Unusually for a polydisperse system, the critical point is found to lie at or very close to the extremum 
of the coexistence region in all cases. We outline an argument showing that such behaviour will occur 
whenever size polydispersity affects only the range, rather than the strength of the inter-particle 
interactions. 
PACS numbers; 64.70Fx, 68.35.Rh 
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I. INTRODUCTION AND BACKGROUND 

A fluid is termed polydisperse when its constituent par- 
ticles are not all identical, but exhibit continuous vari- 
ation in some physical attribute (a, say) such as size, 
shape, charge etc. The state of the system is then quan- 
tifiable in terms of a density distribution p{(j) , measuring 
the number density of particles of each a [1, 2]. Conven- 
tionally, one identifies two distinct classes of polydisper- 
sity: variable and fixed. Variable polydispersity pertains 
to systems such as ionic micelles, oil- water emulsions and 
blood [3-7] where the degree of polydispersity (as mea- 
sured by the form of p(cr)) can change under the influ- 
ence of external factors such as pressure, temperature or 
pH. By contrast, for systems exhibiting fixed polydisper- 
sity, the relative proportions of particles of different a are 
prescribed by the process of their chemical manufacture. 
Examples in this category include colloidal dispersions, 
liquid crystals and polymer solutions [8-10]. For such 
systems, the shape of p{a) is constant, only its scale can 
vary depending on the quantity of solvent present. 

In the present work we focus our attention on fiuids ex- 
hibiting fixed polydispersity. A central issue in such sys- 
tems is the nature of their phase behaviour. It has long 
been appreciated that this can be considerably more com- 
plex than that of monodisperse systems [11]. Indeed, in 
some cases the presence of polydispersity is predicted to 
engender completely new phases, such as the cascades of 
demixing transitions that occur in length-disperse poly- 
mers [12] or the splitting of the triple point in polydis- 



perse liquid-crystal polymers [13]. But even when poly- 
dispersity is not a prerequisite for the existence of a phase 
transition, (e.g. colloid-solvent demixing in colloids), it 
can considerably enrich the character of the transition. 
The reason for this is fractionation: at phase coexistence, 
particles of each a may partition themselves unevenly be- 
tween two (or more) coexisting 'daughter' phases as long 
as-due to particle conservation-the overall composition 
p^{o) of the 'parent' phase is maintained. The conse- 
quences of fractionation for phase diagrams can be dra- 
matic. For example, the conventional liquid-gas binodal 
of a monodisperse system (which connects the ends of 
tie-lines in a density-temperature diagram) splits into a 
'cloud' and a 'shadow' curve. These give, respectively, 
the density at which phase coexistence first occurs and 
the density of the incipient phase. The curves do not 
coincide because the shadow phase in general differs in 
composition from the parent. Further fractionation ef- 
fects are evident when the phase diagram is represented 
in terms of external fields, such as temperature and pres- 
sure. In contrast to monodisperse systems, for which co- 
existence occurs along a line in the pressure-temperature 
projection of the phase diagram, in the presence of poly- 
dispersity, this line is broadened into a coexistence region 
[14, 15]. 

Only recently has experimental work started to clar- 
ify in a systematic fashion the generic consequences of 
polydispersity for phase coexistence properties [16-18], 
and fundamental questions remain unresolved. Consider, 
for instance, a prototype model comprising a system of 



spherical particles exhibiting a short ranged repulsion 
and a long ranged attraction. Even in such a basic sys- 
tem, polydispersity may be manifest in different forms. 
The simplest case is that of a purely size-disperse fluid in 
which the sole variation amongst particles occurs in their 
diameters. Alternatively the amplitude of interparticle 
interactions might vary between particles. More realis- 
tically, one should almost certainly consider a coupled 
combination of both size and amplitude polydispersity 
i.e. the strength with which a particle interacts with a 
neighbor depends on the sizes of both. To date, how- 
ever, there is relatively little known about the individual 
effects of size and amplitude variations in systems ex- 
hibiting fixed polydispersity in terms of their influence on 
phase behaviour. For instance it is unclear just how the 
liquid- vapor critical temperature and density are affected 
when polydispersity is introduced into a model; how the 
critical point shifts depend on the degree and functional 
form of the polydispersity; what is the extent and na- 
ture of fractionation effects in the coexistence region and 
what is the degree of associated coexistence curve broad- 
ening. In view of this, there clearly exists a need for 
systematic computational and theoretical investigations 
of liquid-vapor phase behaviour in model polydisperse 
fluids. In the present paper, we address this need by 
combining Monte Carlo (MC) simulation and analytical 
calculations to investigate liquid-gas phase equilibria in 
size-disperse LJ fluids 

As regards previous studies of phase equilibria in poly- 
disperse system, we know of no prior simulation studies 
of fluids subject to flxed polydispersity (although several 
have been performed for the variable case [19-22]). Some 
analytical studies do exist, however, and these typically 
seek to calculate the system free energy as a function of a 
set of density variables. Unfortunately, this task is com- 
plicated by the fact that the requisite free energy /[p(o-)] 
is a functional of p{a), and therefore occupies an inflnite 
dimensional space. Consequently the standard tangent 
plane approach for identifying phase boundaries becomes 
excessively unwieldy, both conceptually and numerically, 
and one is normally obliged to resort to approximation 
schemes [11]. For sufficiently simple model free energies 
which generalize the van der Waals (vdW) approximation 
to polydisperse systems, a direct attack on the solution 
of the phase equilibrium conditions is nevertheless some- 
times possible, see e.g. [23, 24]. The reason for this is that 
such models are normally "truncatable" so that the phase 
equilibrium conditions can be reduced to nonlinear equa- 
tions for a finite number of variables (see sec. IV below). 
This approach has been applied to the study of phase 
separation in fluids exhibiting separate size and ampli- 
tude polydispersity, yielding predictions for the cloud and 
shadow curves and critical parameters as a function of 
polydispersity [15]. 

An approach which more systematically exploits the 
advantages of truncatable models is the moment free en- 
ergy method [11, 12]. This approximates the full free 
energy appropriately in terms of a "moment free energy" 



which depends on a small number of density variables, 
thereby permitting the use of standard tangent plane 
construction to locate phase boundaries. It thus restores 
much of the standard geometrical insight to the phase 
equilibrium analysis, while at the same time being com- 
putationally efficient. In particular, the method delivers 
(for the given free energy) exact results for the location 
of critical points and the cloud and shadow curves; other 
aspects of phase coexistence can be calculated within a 
systematically refinable approximation scheme that even- 
tually allows one to locate the exact solution of the phase 
equilibrium conditions. The moment free energy has 
been applied to the study of phase behaviour in systems 
ranging from polydisperse hard rods to the freezing of 
hard spheres [11, 25], and we shall deploy it again in the 
present study. 

The paper is arranged as follows [26]. We begin in 
sec. II by defining the form of the parent density distri- 
bution that we have elected to study. Next, in sec. Ill, 
we introduce our simulation model and outline the bat- 
tery of specialized MC techniques necessary to cope with 
fixed polydispersity and phase coexistence within a grand 
canonical ensemble framework. The formalism underpin- 
ning the MFE method, together with the form of the 
model free energy to which it has been applied, are de- 
scribed in sec. IV. Turning then to our results, sec. V 
compares the phase behaviour emerging from the MFE 
calculations with that determined by the simulations. Fi- 
nally, sec. VI presents a discussion of our findings and 
highlights some issues meriting further investigation. 



II. PARENT DISTRIBUTION 

We consider systems characterized by a single continu- 
ous polydisperse attribute a, with associated density dis- 
tribution p{(7). Under conditions of fixed polydispersity, 
the shape of p{a) is imposed. However, its scale can vary 
depending on the amount of solvent present. Following 
convention, we define a "parent" distribution p'^{<y) as 
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where uq — N/V is the overall particle number density, 
while f{cr) is a nominated normalized shape distribution. 
Since /^''(cr) may vary only in terms of its scale no, the 
system is constrained to traverse a so-called dilution line 
in the full phase space of possible compositions. The 
value of no thus serves to parameterize this line. 

The term "parent" is adopted to emphasize the fact 
that although (for a given point on the dilution line) the 
density distribution is fixed globally, the system might 
nevertheless phase separate into coexisting "daughter" 
phases whose density distributions differ from the parent. 
The daughter phase distributions arc then related to the 
parent via a simple volumetric average: 



x(i)p(i)(a) 



-a;(2)p(2)(^) 



.°(a) 



(2) 



where Vx^^'^ and Vx^"^"^ are the volumes of phases (1) and 
(2) respectively, with V the sample volume. Note that on 
the cloud curve marking the onset of phase coexistence, 
the volume fraction of the incipient (shadow) phase is 
negligible and thus the density distribution of the cloud 
phase coincides with that of the parent. 

In this work, we will principally be concerned with 
unimodal shape distributions of the Schulz form: 



fW) = -r 
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a exp 
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(3) 



Here tr = 1 is the average particle diameter, while z is a 
parameter controlling the width of the distribution and 
hence the dimensionless degree of polydispersity. The 
latter is conventionally defined as the standard deviation 
of p(cr), normalized by the mean: 
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III. SIMULATION ASPECTS 

Possibly the simplest MC simulation strategy for ob- 
taining the thermodynamic properties of a polydisperse 
fluid is to populate a simulation box of fixed volume V 
with an appointed number of particles N whose sizes are 
drawn from the desired p'^{cr) [27]. Operationally, how- 
ever, use of such a canonical ensemble is far from optimal 
because it only samples a single realization of the possi- 
ble configurations of particle sizes. Moreover the fixed 
overall particle number prevents effective study of phase 
separation phenomena. 

Experience with the simulation of monodisperse flu- 
ids has shown that use of the grand canonical ensemble 
(GCE) is highly effective for studying fluid phase tran- 
sitions [28-30]. Its application in the context of poly- 
disperse fluids retains the benefits of the monodisperse 
case. Moreover, it provides the key to improved sam- 
pling of the density distribution /o(cr). This is because 
p{a) can fiuctuate as a whole, thus capturing many dif- 
ferent individual realization of the ensembles of particle 
sizes and hence catering naturally for fractionation ef- 
fects. Notwithstanding these advantages, however, the 
GCE might appear (at first sight) unsuitable for the pur- 
pose of traversing a dilution line. This is because the en- 
semble averaged density distribution p{a) ostensibly lies 
out-with the direct control of the simulator, its form be- 
ing instead determined by the imposed chemical potential 
distribution /i(cr). Nevertheless, as we have shown pre- 
viously [31], it turns out to be possible to adapt /x(cr) in 
such a way as to realize a specific desired form of p(o') 
for any temperature of interest. By complementing this 
approach with extended sampling MC techniques, phase 
coexistence properties can be studied along a dilution 



line. The remainder of this section describes the neces- 
sary techniques and the strategy for their implementa- 
tion. 



A. Model and algorithm 

We consider a fluid of particles interacting via a pair- 
wise potential of the Lennard- Jones (LJ) form: 



U{r,. 



4e 



(5) 



with a-ij = ((Ti + (Tj)/2. A cutoff was applied to the 
potential for particle separations r^ > 2.5<Tij. The poly- 
dispersity enters solely through the continuous distribu- 
tion of diameters Ci, which were assigned the Schulz form 
Eq. (3). 

The GCE Monte Carlo algorithm employed to study 
this model invokes four types of operation: particle dis- 
placements, deletions, insertions, and resizing. The par- 
ticle diameter a is treated as a continuous variable in 
the range < ct < CTc, with CTc a cutoff which we choose 
as CTc = 3 unless otherwise specified. However, distri- 
butions defined on a such as the instantaneous density 
p{<t), and the chemical potential h{(t), are represented as 
histograms defined by discretising the permitted range 
of a into 120 bins. For further details concerning the 
simulation algorithm, as well as the structure, storage 
and acquisition of data, we refer the interested reader to 
rcf. [31]. 

The primary observables with which we shall be con- 
cerned are the instantaneous density distribution p{a) 
and its ensemble averaged form p{a). We also accumu- 
late the distribution of several quantitites which do not 
depend explicitly on a. These are the fluctuating over- 
all number density p{n), where n — J dap{a), the dis- 
tribution of the fluctuating volume fraction p{ri), with 
rj — (tt/6) J daa^ p{a), and the distributions of the con- 
flgurational energy p{E) and the energy per particle p{u) 
with u = E/N. 

For simulations of fixed polydispersity at some given 
temperature T, one requires that form of the chemi- 
cal potential distribution p{a) for which p{a) matches 
some target, namely the prescribed parent p^{<y)- The 
task of determining the desired p.{(T) is complicated by 
the fact that it is an unknown Junctional of the parent 
(i.e. p{<y) = /i[p°((T)]). Recently, however, a computa- 
tional tool has been developed that efficiently overcomes 
this difficulty. The non-equilibrium potential refinement 
(NEPR) scheme [32] enables the efficient iterative de- 
termination of p,[p^(a),T], from a single simulation run, 
and without the need for an initial guess of its form. To 
achieve this, p,{(t) is continually updated (in the course 
of a simulation run) in such a way as to minimize the 
deviation of the instantaneous density distribution p{a) 
from the desired parent form. This procedure realizes 
a non-equilibrium steady state for which p{(j) — p^((y)- 



However, since tuning ii{a) 'on-the-fly' in this manner 
violates detailed balance, the form of fJ,{a) thus obtained 
is not the equilibrium solution one actually seeks. Nev- 
ertheless by performing a series of iterations in which 
the degree of modification made to /i(cr) at each step is 
successively reduced, one can drive the system towards 
equilibrium, thereby obtaining the correct fi[p'^{(7),T]. 

Notwithstanding the utility of the NEPR method for 
solving the inverse problem ij,[p^{(7)], its deployment to 
determine the form oi ii{a) all along a dilution line would 
represent a considerable computational endeavor. Fortu- 
nately, however, this is not necessary. A single appli- 
cation of the NEPR method, to obtain p,{a) at a state 
point on the dilution line close to the region of inter- 
est, suffices to bootstrap an efficient tracking procedure 
based on Histogram Extrapolation [33] (HE) techniques. 
Use of HE permits a stepwise exploration of the dilution 
line to lower or higher densities. The essential idea is 
to statistically reweight the data for p{a) measured at 
some dilution line point (given by the value of no) in or- 
der to estimate the form of p{a) corresponding to some 
other not-too-distant value of no — n^. This is achiev- 
able, within the HE scheme, by minimizing a cost func- 
tion measuring the deviation of the extrapolated form 
of p{a)[iJ,{a)] from a target corresponding to the desired 
parent n'^f{(j). Full details of the tracking procedure can 
be found in ref. [31]. 



B. 



Coexistence curve mapping strategy 



Simulation studies of phase coexistence present distinc- 
tive computational challenges. Principal among these is 
the large free energy (surface tension) barrier separat- 
ing the coexisting phases. In order to accurately locate 
coexistence points, a sampling scheme must be utilized 
which enables this barrier to be surmounted [34]. One 
such scheme is multicanonical preweighting [35], which 
utilizes a weight function in the MC acceptance proba- 
bilities, in order to encourage the simulation to sample 
the interfacial configurations of low probability. At a 
given coexistence state point, the requisite weight func- 
tion takes the form of an approximation to the inverse of 
the distribution of the fluctuating number density, p{n) . 
While specialized techniques allow determination oip{n) 
from scratch, in situations where one wishes to track 
a fluid-fluid phase boundary prior determination of a 
weight function is unnecessary, provided one commences 
from the vicinity of the critical point where the barrier 
to inter-phase crossings is small [30]. Data accumulated 
here can be used (together with HE) to provide estimates 
of suitable multicanonical weight functions at lower tem- 
peratures [28] where the barrier height is greater. 

We seek the intersection of the dilution line with the 
temperature dependent coexistence region. The latter 
is delineated by the cloud curves which herald the on- 
set of coexistence when approached from the respective 
pure phases. The cloud curves (and their corresponding 



shadow curves) were obtained as follows. The dilution 
line tracking procedure was bootstrapped at a gas phase 
state point on the dilution line at a moderately low tem- 
perature by using the NEPR method [32] to determine 
the form of p{a). Starting from this point, the dilution 
line was then followed towards increasing density (with 
the aid of HE) until the gas spontaneously liquefied. Hav- 
ing estimated the location of a coexistence state point 
in this manner, the temperature was increased in steps 
(whilst remaining on the dilution line) until the density 
difference between the gas and the spontaneously formed 
liquid vanished, signalling the proximity of the critical 
point. Finite-size scaling methods [28] were then used to 
furnish more precise estimates for the critical parameters. 

Having located the critical point, a detailed mapping 
of the cloud and shadow curves was performed. The key 
to achieving this is the form of the fluctuating overall 
number density, p{n). The gas phase cloud point (incipi- 
ent liquid phase) corresponds to the situation where p{n) 
is bimodal, but with vanishingly small weight in the liq- 
uid peak. Under these conditions, the position of the low 
density gas peak provides an estimate of the gas phase 
cloud density, while that of the liquid peak gives the liq- 
uid shadow density. The converse is true for the liquid 
phase cloud point and its gas shadow. Determining the 
cloud and shadow points as a function of temperature 
yields the cloud and shadow curves. One tracks the gas 
and liquid cloud curves (and their shadows) in a step- 
wise fashion downwards in temperature from the critical 
point, using HE to negotiate each temperature step. HE 
yields estimates not only for the form of ii{a) on the 
cloud curve at the next temperature, but also the requi- 
site multicanonical weight function. It should be pointed 
out that while the positions of the peaks in p{n) provide 
an accurate estimate of cloud and shadow points at low 
temperatures, this breaks down near the critical point 
due to flnite-size effects [28]. Thus a naive extrapolation 
of our curves to their intersection point will tend to over- 
estimate the critical temperature. However, our indepen- 
dent determination of the critical point using flnite-size 
scaling methods is considerably more accurate. 

In order to explore the coexistence region one must 
traverse the dilution line connecting the cloud points at 
some desired subcritical temperature. The strategy for 
doing so commences at one cloud point and entails con- 
structing a series of simulations in which the value of no 
is successively increased. The requisite multicanonical 
weight function for each simulation state point is again 
found with the help of HE from the previous state point. 



IV. MOMENT FREE ENERGY METHOD 

A. Model 

To describe our system theoretically, we need a suitable 
free energy; this is a functional of the density distribution 
p{a) and an ordinary function of temperature. The free 



energy (density) of a polydisperse system can generally 
be decomposed as 



/= dap{a)[\npia)-l] + r 



(6) 



where the first part is the free energy of an ideal polydis- 
perse mixture. To find a suitable model the excess free 
energy f^^, we approximate the repulsive part of our LJ 
interaction as completely hard. For the corresponding 
contribution to f°^ we use the BMCSL approximation 

/bmcsl [^^ 3*^' ^'''l- ^^ ^^"^ monodisperse case this re- 
produces the Carnahan-Starling equation of state [38]; 
in the general polydisperse case it is a function of the 
moments up to third order of the density distribution, 
Pi = Jdaa^p{a) [i = 0...3). We then treat the at- 
tractive part of the interaction potential in the simplest 
possible way, by adding a quadratic van der Waals term 
to f^^. Using the fact that the interaction volume of two 
particles of diameters a and a' scales as [a + cr'Y , this 
gives for the overall excess free energy 



important feature of our excess free energy is that it is 
truncatable, i.e. depends only on a finite set of moments 
Pi — Jda Wi(cr)p{a) defined by weight functions Wi(cr); 
in our case Wi{a) = a^. In coexisting phases, the chem- 
ical potentials fj,{a) and pressure P must be equal. The 
former are, by differentiation of Eq. (6), 

m('^) = Sz^ = In p(a) + ^ ^,rw. (a) , ^^f - ^^" 
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(9) 
while the pressure is given by the Gibbs-Duhem relation 



P = -f + J da fi{a)p{a) = p, - /- + ^ 
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To the conditions of equality of chemical potentials and 
pressure we need to add the requirement of conservation 
of particle number for each species ct, which reads 
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,r = /^MCSL-^ ldad<j'p{a)p{<j'){a + a'f (7) 
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where t is an appropriate dimensionless temperature. 
We note in this context that for our theoretical calcu- 
lations we measure all diameters in units of ct, all densi- 
ties in units of the inverse volume of a reference particle, 
l/[(7r/6)(T^] and all energies in units of fceT. The free 
energy densities in Eq. (6) and Eq. (8), for example, are 
in units of fcBT/[(7r/6)(7^]. Also, p^ is just the volume 
fraction of particles, while po is a scaled number density, 
Po = (7r/6)(T^n. With these conventions, the BMCSL 
contribution to the excess free energy takes the form 
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Given the approximate character of our model free en- 
ergy, it would not make sense to try to scale the tem- 
perature parameter t precisely to the temperature in the 
simulations. Our main aim is to study whether our f^^ 
can reproduce the qualitative trends observed in the sim- 
ulations. While still somewhat crude, our /"'^ is better 
suited to this task than previous versions [15] because 
it incorporates polydispersity not only into the attrac- 
tive contribution, but also into the hard core reference 
system. 



B. Moment free energy 

Our computational approach for determining the phase 
behaviour of the model defined by Eq. (8) is based on the 
moment free energy (MFE) method. We give a brief out- 
line here; details can be found in [11, 12, 39, 40]. The 



where a = 1 , . . . , p labels the phases (compare equa- 
tion Eq. (2)). One then finds from equality of the p{cr), 
Eq. (9), together with particle conservation Eq. (11) that 
the density distributions in coexisting phases can be writ- 
ten as 
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Here the A^ must obey 



Af^=-Mf^"'^ + c. 



(13) 



and the Ci are undetermined constants that do not affect 
the density distributions Eq. (12). One can fix them e.g. 
by requiring all the Aj" in one of the phases to be zero. 
A little reflection then shows that Eq. (13) together with 
Ea ^^"^ ~ 1 ^'^'^ ^^^ equality of the pressures Eq. (10) 
in all phases give a closed system of nonlinear equations 
for the p{M + 1) variables A^" and x^"-\ A solution 
can thus, in principle, be found by a standard algorithm 
such as Newton-Raphson. Generating an initial point 
from which such an algorithm will converge, however, is 
still a nontrivial problem, especially when more than two 
phases coexist and/or many moments pi are involved. 
Furthermore, the nonlinear phase equilibrium equations 
permit no simple geometrical interpretation or qualita- 
tive insight akin to the construction of phase diagrams 
from the free energy surface of a finite mixture. 

The moment free energy addresses these two disadvan- 
tages. To construct it, one starts by modifying the free 
energy decomposition Eq. (6) to 
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+ ri{p^}) (14) 



In the first (ideal) term, a normalizing factor i?(cr) has 
been included inside the logarithm. This has no effect on 
the exact thermodynamics because it contributes only 
terms linear in /o(cr), but will play a central role below. 
One can now argue that the most important moments 
to treat correctly in the calculation of phase equilibria 
are those that actually appear in the excess free energy 
f^^iiPi})- Accordingly, one imposes particle conserva- 
tion Eq. (11) only for the pi, but allows it to be violated 
in other details of the density distribution p(a) which do 
not affect the pi . These "transverse" degrees of freedom 
are instead chosen to minimize the free energy Eq. (14), 
and more precisely its ideal part since the excess con- 
tribution is a constant for fixed values of the pi. This 
minimization gives 



p{a) = R{(7) exp 
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where the Lagrange multipliers Xi are chosen to give the 
desired values of the moments 



Pi = da Wi (a) R{a) exp 
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The corresponding minimum value of / as given in 
Eq. (14) then defines the moment free energy (MFE) 



fmomilPi}) = i^XiPi - 
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Since the Lagrange multipliers are (at least implicitly) 
functions of the moments, the MFE depends only on 
the Pi. These can now be viewed as densities of "quasi- 
species" of particles, allowing for example the calculation 
of moment chemical potentials [12] 
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and the corresponding pressure P — ^^ piPi — /mom 
which turns out to be identical to the exact expression 
Eq. (10). A finite-dimensional phase diagram can thus 
be constructed from /mom according to the usual tan- 
gency plane rules, ignoring the underlying polydisperse 
nature of the system. Obviously, though, the results now 
depend on R{a). To understand its influence, one notes 
that the MFE is simply the free energy of phases in which 
the density distributions p{a) are of the form Eq. (15). 
To ensure that the parent phase is contained in the fam- 
ily, one normally chooses its density distribution as the 
prior, R{a) = p^{cr); the MFE procedure will then be 
exactly valid whenever the density distributions actually 
arising in the various coexisting phases are members of 
the corresponding family 



pM =P°(CT)exp 
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It is easy to show from Eq. (12) that this condition holds 
whenever all but one of a set of coexisting phases are 
of infinitesimal volume compared to the majority phase. 
Accordingly, the MFE yields exactly the onset of phase 
of coexistence, conventionally represented via cloud and 
shadow curves as discussed above. Similarly, one can 
show that spinodals and critical points are found ex- 
actly [12]. 

For coexistence involving finite amounts of different 
phases the MFE only gives approximate results, since 
different density distributions from the family Eq. (19), 
corresponding to two (or more) phases arising from the 
same parent p" (a) , do not in general add to recover the 
parent distribution itself. Moreover, from Gibbs' phase 
rule, a MFE depending on M moments will not predict 
more than M + 1 coexisting phases, while we know that 
a polydisperse system can in principle separate into an 
arbitrary number of phases. Both of these shortcomings 
can be overcome by including extra moments within the 
MFE; this systematically increases the accuracy of any 
calculated phase splits [12]. By choosing the weight func- 
tions of the extra moments adaptively, the properties of 
the coexisting phases can then be predicted with in prin- 
ciple arbitrary accuracy [12, 41]. Importantly for us, the 
results can in fact be used as initial points from which a 
solution of the exact phase equilibrium problem can be 
converged successfully [42, 43]. This is the technique that 
we use here to obtain results in the coexistence region. 



V. RESULTS 

Here we report the results of applying the simulation 
and MFE methods described in sec. Ill and IV to obtain 
the phase behaviour of size disperse fiuids. Our findings 
are separated into two parts, those for the sub critical re- 
gion (sec. V A) and those for the critical region (sees. VB 
and VC). With regard to the representation of results, 
the MFE calculations are essentially exact for the given 
model free energy; symbols indicate the temperatures at 
which individual calculations were performed and lines 
are merely guides to the eye. For the simulation results, 
unless otherwise indicated, statistical errors do not ex- 
ceed the symbol sizes. Again, lines are guides to the eye. 



A. The subcritical region 

Computational complexity limits our investigation of 
liquid-gas coexistence to the case of a fluid described by 
a size distribution of the Schulz form, Eq. (3) having 
width parameter z = 5. This choice of z corresponds to 
a rather large {S = 40.7%) degree of polydispersity. In 
both the simulations and the MFE calculation, f{a) was 
truncated at CTc = 3.0 and renormalized appropriately. 
The simulations were all performed using a simulation 
box of volume V = 11390(7^. In what follows we quote 



all simulation quantities in the standard dimensionless 
LJ units (e.g. t = kBT/Ae). 
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FIG. 1: (a) Cloud and shadow curves in the n — T plane, 
as obtained from the simulations, (b) The corresponding 
prediction of the MFE calculations for a van der Waals model. 
Also shown for comparison in each case is a portion of the 
coexistence binodal for the monodisperse limit. 

We commence by presenting (in fig. 1) a comparison 
of the simulation measurements of the cloud and shadow 
curves in the n ~ t projection, with those of the MFE 
calculations. Apparent in this representation is a stark 
separation of the cloud and shadow curves. Furthermore, 
the whole phase diagram is considerably shifted with re- 
spect to that of the monodisperse fluid (itself determined 
previously in the case of fig. 1(a) in ref. [28]). Specifically, 
one observes that the critical point occurs at a consider- 
ably higher temperature than in the monodisperse limit. 
We note that this particular finding contrasts with that 
of a previous theoretical study of a size-disperse van der 
Waals fluid [15], which predicts a suppression of the crit- 
ical temperature with respect to the monodisperse limit. 

The order of cloud and shadow curves commonly ob- 
served with increasing density in many polydisperse sys- 
tems (see e.g. [44]) is cloud-shadow-cloud-shadow. By 
contrast, the order shown in flg. 1 is cloud-shadow- 
shadow-cloud. Interestingly, however, the order reverts 
to the standard pattern if one plots the data in terms of 
the volume fraction 77 = (ir/d) J daa^p{a), rather than 
the overall number density, as shown in flg. 2. Moreover, 
one sees that in the rj—t representation the differences be- 
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FIG. 2: (a) The simulation, and (b) MFE data of fig. 1, 
both re-expressed in terms of the volume fraction rj = 

{tv/6) J da a'* p{a). 



tween cloud and shadow phase properties become much 
less pronounced. This is particularly evident in the MFE 
results, flg. 2(b), for which the cloud and shadow curves 
almost coincide, though a close inspection reveals that 
they are in fact distinct and occur in the same order as 
observed in the simulations. Similar flndings pertain to 
the energy per particle measured in the cloud and shadow 
phases (flg. 3), though the qualitative agreement in the 
shape of the curves between simulations and theory is 
less good here. One reason for this that in our model 
free energy the repulsive interactions, which are mod- 
elled as hard, do not contribute to the energy u. In the 
MFE calculations we therefore determine u as the value 
of the attractive van der Waals term from (Eq. 8). In 
the simulations, on the other hand, both attractive and 
repulsive interactions contribute to the measured value 
of u. 

With regard to the critical point parameters, we note 
that while the critical number density of the polydisperse 
fluid is considerably less than its value in the monodis- 
perse limit, the simulation estimates of the critical vol- 
ume fraction for the mono- and polydisperse fluid agree 
to within error. Moreover, in both the simulations and 
the MFE results, the critical point is located at or at 
least extremely close to the maximum of the cloud and 
shadow curves. We see no evidence for portions of the 
cloud or shadow curves lying above the critical temper- 
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FIG. 3; The energy per particle on the cloud and shadow 
curves as a function of temperature as obtained from (a) the 
simulations and (b) the MFE calculation. 



ature, which would generically be expected for polydis- 
perse systems [11]. We postpone further discussion of 
these observations to sec. VC. 

Not all distinctions between cloud and shadow curves 
can be disguised by simply recasting the data in terms 
of ?7, rather than n. At temperatures significantly below 
criticality, we observe considerable broadening of the co- 
existence curve in the space of /.t(cr). This is evident in 
fig. 4 which shows the form of ii{a) at the respective cloud 
points (marking the boundaries of the two phase region) 
for the lowest temperature studied in the simulations, 
t = 0.91tc- Such broadening does not occur in monodis- 
perse systems — coexistence occurs at a single value of 
the chemical potential, not a range of values. The effect 
is surprisingly large, notwithstanding the high degree of 
polydispersity of the parent. Indeed, in simulation terms 
the respective cloud points are so far separated in phase 
space that to connect them directly (via a route crossing 
the phase boundary) required some twenty overlapping 
simulations — three times the number required to connect 
the cloud point to the critical point at this temperature. 

Turning to the form of the density distributions in 
the shadow phases, an example of the scale of the dif- 
ference between these and the parent is given in fig 5 
for t — 0.91ic- The data show that at the gas phase 
cloud point, larger particles preferentially occupy the liq- 
uid shadow phase. Conversely at the liquid phase cloud 
point, there is a predominance of smaller particles in the 
gas shadow phase. Clearly the scale of these fraction- 
ation effects is significant. We note that the principal 
fractionation effect is a change in the mean particle di- 
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FIG. 4: The form of the chemical potential distribution ^{o) 
at the gas phase and liquid phase cloud points at f = 0.91tc. 
(a) Simulation, (b) MFE calculations. 



amcter. The simulations give for the mean diameter in 
the liquid shadow [45] {a) — 1.167(3), while in the gas 
shadow it is {a) = 0.863, both to be compared to that 
of the parent, {a) — a — I. The changes in the shape of 
the distribution (as indicated by polydispersities of the 
shadow phases), are rather smaller. For this temperature 
(t ~ O.Qltc), we find in the simulations S ss 43% for the 
liquid shadow, and S w 38.5% for the gas shadow, to be 
compared with a parent polydispersity oi 6 — 40.7%. 
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— gas shadow 
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FIG. 5: (a) Simulation estimates of the normalized daugh- 
ter phase density distributions on the shadow curves at t = 
0.91tc. Also shown for comparison is the parent distribution. 
(b) The corresponding results of the MFE calculations. 



The temperature dependence of the mean particle size 
and degree of polydispersity S are shown in fig. 6. ft 
is interesting to note that while the value of S for the 
liquid shadow initially increases strongly with decreasing 
temperature, on further reducing temperature it subse- 
quently bends back to lower values of S. (Analogous non- 
monotonic behaviour is observed in the gas shadow.) Our 
data (in conjunction with Eq. (4)) show that the origin of 
this effect lies in the fact that while (a) increases mono- 
tonically with decreasing temperature, the standard devi- 
ation of the size distribution in the liquid shadow phase 
first increases strongly as t is lowered from its critical 
value, but then saturates to a constant value for temper- 
atures below about 0.9<c- 
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in the liquid phase density is considerably greater than 
that for the gas phase. 
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FIG. 7: Distribution of the fluctuating number density p{n) 
at a selection of values of the parent density no between the 
cloud points at f = 0.91ic. 

Finally in this section, we plot in fig. 9 the measured 
form oi p{n) for state points spanning the coexistence re- 
gion, expressed on a logarithmic scale. This representa- 
tion exposes the magnitude of the surface tension barrier 
separating the coexisting phases, which is proportional 
to the logarithm of the peak to trough ratio of the re- 
spective distributions [46] . One observes that the surface 
tension of the liquid shadow at the gas phase cloud point 
is considerably less than that of the gas shadow at the liq- 
uid phase cloud point. This finding reflects the fact that 
the liquid shadow comprises particles that are on average 
larger than those of both the parent (cf. fig. 5) and the 
gas shadow. It follows that the number of particles per 
unit surface area, and hence the free energy cost of the 
interface, is smaller between the gas cloud and the liq- 
uid shadow than between the gas shadow and the liquid 
cloud. 



FIG. 6: Simulation estimates of the degree of polydispersity 
(a) and the mean particle diameter (b) on the shadow curves. 
(c,d) The corresponding results of the MFE calculations. 

In order to explore the coexistence region separating 
the cloud curves, we have scanned the dilution line at 
t — 0.91tc- The simulations yield the distribution of 
the fluctuating number density p{n), the form of which 
is shown in fig. 7 for a selection of values of hq span- 
ning the coexistence region. One observes that in addi- 
tion to a transfer of weight between the gas and liquid 
phase peaks as hq is increased from its value at the gas 
phase cloud point, both the densities of the gas and liq- 
uid peaks shift markedly. This shift of course reflects 
the non-coincidence of the cloud and shadow densities 
for a given phase. We note further (fig. 8) that the den- 
sity shifts are non-linear in ng and that the overall shift 



B. Critical point shifts 

Here we enquire how the critical point parameters of 
the size disperse fluid depend on the degree of polydis- 
persity S. To this end we have employed simulation to 
investigate four Schulz distribution with width param- 
eters (cf. Eq. (3)) z — 50, 25, f and 5, corresponding 
to polydispersities S = 14%, 19.6%, 30.1%, 40.7% respec- 
tively. For z — 5 and z = 10, the size distribution was 
truncated at CTc ~ 3.0, while for z = 25 and z = 50 a 
cutoff ac = 2.0 was used. 

Standard flnite-size scaling (FSS) methods were used 
to estimate the location of the critical point in each case 
studied. The FSS methodology exploits the critical point 
scale invariance of the fluctuation spectrum of quanti- 
ties such as the number density, as expressed through 
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FIG. 8: Dependence of the gas and liquid densities on the 
overall parent number density no. (a) Simulation results ex- 
tracted from fig. 7. (b) The corresponding results from the 
MFE calculations at t = O.Qlic- 




FIG. 9: The data of fig. 7 plotted on a logarithmic scale in 
order to expose the scale of the free energy barrier separating 
the pure phase gas and liquid states in the coexistence region. 



the form of p(n) [28]. For the z = 50 and z — 2h dis- 
tributions, three system volumes were studied, namely 
Y ^ 3375ct3, SOOOct^, 15625a-3 . For the z = 10 and z = 5 
distributions, somewhat larger simulation cells were re- 
quired, the values chosen being V = 113905-'^, 27000ct'^, 
and V = 52734ct'^. We note that the largest of these 
systems contained over 10'' particles at criticality. 

In order to compare with the simulation results, MFE 
calculations were performed to obtain the critical point 
for width parameters in the range from z = 1 to the 



monodisperse limit z = oo. The resulting comparison 
is shown in fig. 10. We note that while the simulations 
were able to locate the critical temperature tc rather ac- 
curately (the error bar is smaller than the symbol size), 
the estimates of the critical volume fraction rjc carry a 
rather large estimated error because of the difficulty in 
obtaining good statistics for large systems and high de- 
gree of polydispersity. Thus, although we find a strong 
increase in tc with increasing 6, no clear trend is ob- 
served for rjc- The MFE results confirm the increase of 
tc with increasing 6. However, the fractional change is 
substantially less than that observed in the simulations. 
As regards the critical volume fraction, the MFE method 
predicts a significant (5%) decrease in the value of ?7c over 
the range of 6 studied in the simulations. The scale of 
this decrease is, however, larger than any that might be 
considered consistent with the simulation error bars. 



C. Location of the critical point 

An important finding of our investigation of critical 
point shifts is that, for each z studied, the critical point 
occurs at, or very close to, the maximum of the cloud 
and shadow curves. Correspondingly, there is no numer- 
ical evidence of distinct cloud and shadow points at or 
above tc, either in the simulations or the MFE calcu- 
lations. This observation is consistent with earlier re- 
sults obtained from a simpler model free energy [15], but 
clearly calls for an explanation: in a polydisperse sys- 
tem the critical point must be found at an intersection 
of the cloud and shadow curves, and this intersection 
will in general be located below the maxima of the two 
curves [11, 14]. 

One possible way to investigate when the critical point 
is near the top of the cloud (and hence also the shadow) 
curve is to find conditions under which the slope of the 
cloud curve there is zero. This can be done for simple ex- 
cess free energies depending only on a single density mo- 
ment, within the framework of a Landau expansion [14]. 
However, such an approach becomes excessively unwieldy 
when several density moments are involved. We have 
therefore addressed the question in a more indirect way. 
On general grounds the spinodal curve, where a given 
polydisperse phase becomes locally unstable, must lie in- 
side the cloud curve and touch it at the critical point [11]. 
This implies that the critical point is at the maximum of 
the cloud curve exactly when it coincides with the max- 
imum of the spinodal (as happens in monodisperse sys- 
tems) . We can therefore estimate the shift of the critical 
point away from the cloud curve maximum by comparing 
its location with that of the spinodal maximum: as long 
as these two points are close together, the critical point 
should also be close to the cloud curve maximum. 

The advantage of this approach is that both the spin- 
odal curve and the critical point can be calculated rel- 
atively easily within the MFE framework, from local 
properties of the excess free energy surface. Fig. 10 
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FIG. 10: (a) Simulation estimates of the critical tempera- 
ture and volume fraction for Schuiz parent distributions of 
a variety of widths (see text for details) . (b) Corresponding 
predictions of the MFE calculations for Schuiz distributions in 
the range z = 1 to 2 = oo (monodisperse limit). The dashed 
line and empty circle indicate the location of the spinodal 
maximum; see sec. V C. 



shows that the resulting locations of the critical point 
and spinodal maximum arc, though not identical, ex- 
tremely close for the whole range of width parameters 
z studied above. Even for the most polydisperse system 
(z = 1, 5 = l/\/2 K, 70%) the temperature coordinates 
of the two points are indistinguishable on the scale of the 
graph, confirming our earlier observation that the critical 
point is very near the cloud curve maximum. 

To understand under which circumstances such rather 
unusual behaviour can be expected, we have derived 
general expressions for the location of the critical point 
and spinodal maximum in a polydisperse system with a 
smooth (van der Waals-type) excess free energy. Details 
will be given in a forthcoming publication [48]. Within 
an expansion valid for near-monodisperse systems, i.e. 



small polydispersity (5, the leading shifts in the critical 
point and spinodal maximum relative to the monodis- 
perse limit (5^0 turn out to be 0(<5^). The density shifts 
are in general different as expected on general grounds 
for a polydisperse system; see above. (The tempera- 
ture shifts, which are essentially quadratic in the density 
shifts because we are expanding around the maximum 
of the monodisperse spinodal curve, are always identi- 
cal to 0((5^).) However, it turns out that for purely 
size-polydisperse systems, where the excess free energy 
is unchanged if we scale all particle sizes and interparti- 
cle separations by the same factor, the lowest-order den- 
sity shifts of critical point and spinodal maximum exactly 
coincide. Likewise, one can show using the techniques 
of [16, 47] that size polydispersity causes the lowest-order 
shifts in the entire cloud and shadow curves to be iden- 
tical, once they are plotted in terms of volume fraction 
rather than density. This is consistent with our observa- 
tion of the near-coincidence of these curves in fig. 2. 

Thus the perturbative expansion shows that the effects 
which we observe numerically are due to the fact that 
our LJ system has purely size-polydisperse interactions. 
This rationalizes why earlier calculations [15] also found 
the critical point very near the top of the cloud curve, 
even though a much simpler free energy model was used: 
any excess free energy which respects the scale invariance 
of a size-polydisperse system will exhibit such behaviour. 
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FIG. 11: Normalized size distributions in the daughter phases 
on the shadow curves at i « 0.9tc, for a top hat parent dis- 
tribution with polydispersity & = 40%. 

Wc initially thought that the fact that our critical 
points were at the top of the cloud curves was related 
to the size distributions in the shadow phases being es- 
sentially just scaled versions of the cloud (parent) size 
distributions; see fig. 5. Indeed, one can show that for 
any size-polydisperse excess free energy which depends 
only on the density moments po and pi , this scaling holds 
exactly when the parent phase has a Schuiz size distri- 
bution [48]. In such a scenario the cloud and shadow 
curves also coincide exactly in the t] — t representation, 
and the critical point is exactly at their common max- 
imum [48]. However, the same features can occur, to a 
very good approximation, even when there is no scaling 
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FIG. 12: (a,b) Cloud and shadow curves in the n ~ t and 
ri — t planes, for a parent with a top-hat size distribution, as 
obtained by MFE calculations. (c,d) The analogous results 
for a system with added amplitude polydispersity; the critical 
point is now clearly below the maximum of cloud and shadow 
curves. 



link between cloud and shadow size distributions. We 
demonstrate this by considering theoretically a parent 
phase with a top hat size distribution, with the poly- 
dispersity 5 = 40% taken to be the same as our main 
Schulz distribution example. Fig. 11 shows exemplary 
size distributions in the shadow phases: these are now 
clearly different in shape from each other and from the 
parent. Nevertheless (fig. 12(a,b)) the rj — T cloud and 
shadow curves are almost identical and the critical point 
is extremely close to their maximum. 

Finally, the above theoretical arguments suggest that, 
once we move to a case where polydispersity affects not 
only the size but also the amplitude of the interaction 
potentials, the behaviour expected for a generic polydis- 
perse system should be recovered. To confirm this, we 
repeated the MFE calculations for interaction potentials 
Uij, Eq. (5), scaling as the product aiCTj of the sizes of 
the particles involved. This translates into an additional 
factor aa' in the attractive term in the model free en- 
ergy (Eq .7). The results are shown in fig. 12(c,d) for 
the top hat parent: as expected, the critical point is now 
clearly below the cloud curve maximum, and the rj — T 
representations of the cloud and shadow curves no longer 
coincide. Preliminary simulation results for a LJ model 
with combined size and amplitude polydispersity confirm 
this scenario [49]. 



VI. DISCUSSION AND OUTLOOK 

In summary, we have used specialized MC simulation 
and the moment free energy (MFE) method to study 
the liquid-gas phase behaviour of a size-disperse LJ fluid 
having a degree of polydispersity 5 — 40.7%. Cloud and 
shadow curves have been traced and fractionation effects 
quantifled. Surprisingly good qualitative agreement is 
found in almost all aspects of the results between the 
simulations and the MFE calculations for a van der Waals 
model. 

In a related study, we have also obtained the depen- 
dence of the critical point parameters on the degree of 
polydispersity. The simulation show a strong rise in the 
critical temperature with increasing (5, but no clear trend 
in the critical volume fraction rjc- The increase in tc is 
confirmed by the MFE calculation, albeit with a weaker 
magnitude. However, in contrast to the simulation, the 
MFE predicts quite a strong decrease of rjc with increas- 
ing 5. 

Some of our findings for the sub-critical coexistence re- 
gion seemed, at first sight, rather surprising: the cloud 
and shadow curves nearly collapsed onto a common curve 
in the rj — t representation, with the critical point found 
so close to their maximum as to be indistinguishable. 
We outlined a theoretical argument which shows that 
such behaviour is generically to be expected for systems 
with pure size-polydispersity, but not otherwise. This 
was confirmed by, on the one hand, considering a differ- 
ent parent size distribution, with a top hat shape: sim- 
ilar behaviour was observed. The addition of amplitude 
polydispersity, on the other hand, completely changed 
the picture. We demonstrated finally that the close sim- 
ilarity in the shapes of the size distributions of cloud 
and shadow phases is rather peculiar to parents with a 
Schulz size distribution, and not directly linked to the oc- 
currence of the other effects of size polydispersity which 
we observe. 

Given the rather crude nature of the van der Waals 
free energy Eq. (8), used for the MFE calculation, the 
level of qualitative agreement with the simulations in the 
sub-critical region is gratifying. We note, however, that 
previous work treating polydispersity within a van der 
Waals framework [15] did not produce a similar level of 
agreement. We tentatively ascribe this difference to the 
inclusion in the present work of an explicitly polydisperse 
hard-core reference free energy having the well known 
BMCSL form [36, 37]. Use of this reference free energy 
improves the description of packing effects. These are 
expected to be particularly important for the correct de- 
scription of liquid structure in polydisperse fluid because 
of the possibility that small particles can fit into the gaps 
between large ones. We speculate that it is this ability 
to pack more effectively (and the resulting lowering of 
the configurational energy) which is responsible for the 
observed increase in tc with increasing polydispersity, a 
finding which contrasts with the predicted decrease in 
rcf. [15]. 
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One surprise arising from the simulation results is the 
observed magnitude of coexistence curve broadening. In 
simulations of a monodisperse system at coexistence, the 
system fluctuates freely between both phases (assuming 
the provision of an appropriate multicanonical weight 
function) [28]. The form oi p{n) exhibits two peaks, and 
these are visible over a certain range of fi due to finite- 
size smearing of the transition. By contrast the range of 
IJ,{a) (fig. 4) over which gas and liquid phase peaks are 
visible (cf. fig. 7) is 1-2 orders of magnitude greater than 
the finite-size smearing one would expect in a monodis- 
perse system having a similar number of particles. Conse- 
quently a large number of separate, but overlapping (in 
configuration space) simulations are required to bridge 
the coexistence region from pure gas to pure liquid. 

With regard to the general simulation issues raised in 
this study, we have shown that recently developed tech- 
niques (dilution line tracking [31] bootstrapped by the 
non-equilibrium potential refinement method [32]) can 
be successfully employed to map the phase behaviour of 
fluids exhibiting fixed polydispersity. This is achieved 
within what is arguable the most flexible framework for 
the study of fluid phase equilibria, namely the grand 
canonical ensemble. Notwithstanding these advances, 
however, it should be stressed that the simulations them- 
selves are not yet routine: the work reported here con- 
sumed well over 10^ hours of CPU time on a 2 GHz 
PC processor. Such an investment of effort exceeds 
that necessary for obtaining the phase behaviour of a 
monodisperse fluid by a factor of 1-2 orders of magnitude. 



The source of the computational complexity is twofold. 
Firstly, one needs to utilize system sizes whose linear ex- 
tent exceeds a given factor of the largest possible particle 
diameter. Even for moderate polydispersities, this can 
increase the necessary volume by an order of magnitude 
compared to the monodisperse case. Secondly, in order 
to relax the system between sampled configurations, one 
must decorrelate not only the overall density, but also 
the distribution of particle sizes p{a) . The bottleneck for 
the latter process is the largest particles, for which the 
probability of particle insertions and deletion is small, 
and which can thus only be altered by the resizing moves 
whose effect is generally more incremental in magnitude. 
Clearly there is scope for further algorithmic improve- 
ments in this regard. 

As an extension of the present study, work is under way 
to simulate the phase diagram of a Lennard-Jones fluid 
in which the well depth in the LJ interaction potential 
depends on the sizes of the interacting particles. Pre- 
liminary results [49] confirm the prediction of the MFE 
calculations described in sec. VC, i.e. that the cloud 
and shadow curves adopt the standard order and that 
the critical point lies well below the extremum of the 
cloud and shadow curves. 
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